Molecular transport calculations with Wannier functions 
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We present a scheme for calculating coherent electron transport in atomic-scale contacts. The 
method combines a formally exact Green's function formalism with a mean-field description of 
the electronic structure based on the Kohn-Sham scheme of density functional theory. We use an 
accurate plane- wave electronic structure method to calculate the eigenstates which are subsequently 
transformed into a set of localized Wannier functions (WFs). The WFs provide a highly efficient 
basis set which at the same time is well suited for analysis due to the chemical information contained 
in the WFs. The method is applied to a hydrogen molecule in an infinite Ft wire and a benzene- 
dithiol (BDT) molecule between Au(lll) surfaces. We show that the transmission function of BDT 
in a wide energy window around the Fermi level can be completely accounted for by only two 
molecular orbitals. 

PACS numbers: 



I. INTRODUCTION 

The transport of electrons in molecules plays a ma- 
jor role in many scientific areas. In biological systems 
for example electrons have to be transported to or from 
chemically active parts of enzymes. In electrochemistry 
the functioning of an electrochemical cell can be affected 
by the electron transport through molecular layers at the 
electrodes. In molecular electronics the goal is to control 
the electron transport in such detail that electronic com- 
ponents like rectifiers or transistors can be constructed 
from single molecules. Within these different fields a ba- 
sic understanding of many phenomena related to electron 
transport has been obtained, but the challenge still re- 
mains to develop a detailed, quantitative approach to 
accurately determine the transport properties from basic 
quantum mechanical principles. 

In this paper we address this issue in the domain where 
the electronic motion can be regarded as coherent. We 
present a numerical Green's function method for calcu- 
lating the linear response conductance of molecules and 
nano-contacts. The system under investigation is divided 
into a central scattering region and attached leads and 
the electronic structure of both the scattering region and 
the leads is described using Density Functional Theory 
(DFT). To represent the electronic states we use a set 
of localized basis functions consisting of partly occupied 
Wannier functions (WFs). The WFs describe the rele- 
vant states very accurately - in fact they are constructed 
so that they reproduce the results of high-accuracy plane- 
wave based pseudopotential calculations - and at the 
same time they represent a minimal basis set which is 
computationally very efficient. 

In addition to their technical advantages, the Wannier 
functions provide insight into the local chemistry of the 
system and are thus well suited for analysis purposes. 
Being the local analogue of the extended (Bloch) states 
of solid state physics, the Wannier functions formalize 



chemical concepts such as bond types and electron lone 
pairs. The traditional definition of Wannier functions for 
single isolated bands (insulators) or occupied molecular 
orbitals (molecules) does not apply to metallic systems^^. 
In this case an extended scheme involving the inclusion of 
selected unoccupied states through a disentangling pro- 
cedure must be used instead. We present a newly devel- 
oped scheme which solves this problem in a simple and 
very direct way and we show how to link the resulting 
"partly occupied WFs" to the general transport formal- 
ism through the evaluation of the relevant Hamiltonian 
matrices. 

We apply the approach to two different systems. In 
the first we calculate the conductance of a hydrogen 
molecule suspended between two monatomic Pt wires. 
This can be viewed as a simple model of a recent experi- 
ment which investigated the conductance of a molecular 
hydrogen contact between bulk Pt electrodesi. By diago- 
nalizing the Hamiltonian within the subspace spanned by 
the WF of the molecule we obtain renormalized bonding 
and anti-bonding H2 orbitals. The individual contribu- 
tions of these orbitals to the total transmission is then 
studied by directly removing each of the orbitals from 
the basis set. In this way we find that transport is due to 
the anti-bonding state. The second system consists of a 
benzene-dithiol molecule between Au(lll) surfaces. This 
system was among the first single-molecule systems to be 
studied experimentally^ and represents a test system for 
transport calculations. Our results for the conductance 
agree qualitatively with previous calculations, but using 
the analysis technique described above, we furthermore 
show that the transport properties of BDT can be com- 
pletely accounted for by only two molecular orbitals. 

The paper is organized as follows. In Sec. |n] we 
introduce the concept of phase-coherence and give an 
overview of existing methods for transport calculations. 
We then proceed to develop the general Green's func- 
tion formalism underlying the present scheme, and illus- 
trate it through a detailed discussion of the simple case 
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of transport through a single energy level. In Sec. Illll we 
introduce a numerical method to construct partly occu- 
pied Wannier functions and apply it to an isolated ethy- 
lene molecule and an infinite trans-polyacetylene wire. 
In Sec. llVl we connect the Wannier functions to the gen- 
eral transport formalism by showing how to obtain the 
relevant Hamiltonian matrices in a Wannier function ba- 
sis. We put emphasis on issues related to the combina- 
tion of periodic supercells and transport calculations. In 
Sec. |V] we present conductance calculations for a molec- 
ular hydrogen bridge between monatomic Pt wires and 
a benzene-dithiol molecule suspended between Au(lll) 
surfaces. 



II. TRANSPORT THEORY 

We begin this section with a brief discussion of the 
phase-coherent transport regime. In order place our 
method in the large picture we give an overview of previ- 
ously developed transport schemes. The general Green's 
function formalism and the central conductance formula 
are then introduced in the case of non-interacting elec- 
trons moving in a static mean-field potential. Extension 
of the formalism to include ultrasoft pseudopotentials is 
shown to be straightforward. Finally, in order to illus- 
trate the general theory we give a detailed discussion of 
electron transport in the simple but important Newns 
Anderson model. 



A. Phase-coherent electron transport 

Electron transport is said to be phase-coherent if the 
dimensions of the system under study are smaller than 
the phase relaxation length, which is the length over 
which the phase of an electron is randomized due to 
scattering2i4. The phase-relaxation length thus defines 
the length scale on which electron interference can be 
observed. A typical phase relaxation length for Au at 
T = 1 K is around In general, the phase of an elec- 

tron is destroyed by interactions with the dynamic envi- 
ronment such as the other electrons, phonons and mag- 
netic impurities and consequently the phase relaxation 
length increases at lower temperatures. Static scatterers 
such as the potential from the ions in a rigid lattice do 
not destroy the phase since the scattering is elastic and 
the phase change is fixed. 



B. Existing transport schemes - an overview 

During the last decade a variety of numerical methods 
addressing the transport of electrons in atomic-scale sys- 
tems have been developed. Very generally, these methods 
can be divided into wave function and Green's function 
methods. In both cases the description is based on a 



Landauer-Biittiker setup where a central scattering re- 
gion is placed between two ballistic leads which are as- 
sumed to be in thermal equilibrium far from the central 
region. 

In the wave function method one solves directly for the 
scattering wave functions of the system and obtain the 
conductance/current from the elastic transmission coef- 
ficients in accordance with the well known Landauer for- 
mula^. The scattering states can be calculated using 
either the transfer matrix methodi^'^iiliiiiii^j the wave 
function matching methodASi or by solving a Lippman- 
Schwinger equatio n^'^i^^ . These techniques have been 
combined with various approximations for the electronic 
structure of the system, ranging from semi-empirical 
models like extended Hiickel^*^ to fully atomistic descrip- 
tions based on the single-particle Kohn-Sham scheme of 
density functional theory (DFT)2ii^. Most of the wave 
function methods, however, apply an intermediate level 
of approximation where only the electronic structure of 
the central region is resolved in detail while the leads are 
modeled by a free electron gas (jellium)i24iiii2iiiii&. 

As an alternative to the Landauer-Biittiker approach 
the transport problem can be formulated with in the 
non-equilibrium Green's function formalism developed 
by Keldysh^^ and Kadanoff and Baymii. An advan- 
tage of the Green's function (GF) formulation is the pos- 
sibility of extending the description beyond the single- 
particle picture by including e.g. inelastic scattering 
on vibrations or correlation effects through self-energies. 
While the inclusion of such many-body effects is formally 
straightforward, the practical calculation of reliable self- 
energies from first principles is a difficult task which only 
recently has been addressed^*. For non-interacting elec- 
trons, the GF- and Landauer-Biittiker formalisms are 
equivalents^. In this case the main effort of the GF 
method is to calculate the GF of the central region in 
the presence of coupling to the leads. The effect of cou- 
pling to the leads is usually incorporated through non- 
hermitian self energiea^°i^ii'^^.^3,24,25,26,27,28,29,30_ ^he 

various methods mainly differ in the representation of the 
Hamiltonian and self-energy matrices which range from 
tight-binding modelsSSiSi, to fully atomistic DFT de- 
scriptions^iiSiSiiSSiSSiSi. Intermediate descriptions treat- 
ing only the central region at the atomic level and adopt- 
ing a parametrization of the coupling and lead electronic 
structure have also been used22422i2£. 

The main difference among the fully DFT based GF 
methods lies in the basis set used to evaluate the Hamil- 
tonian and thus the Green's function. One limitation 
inherent in the GF formalism is that the basis set should 
consist of fimctions with finite support. As a consequence 
most methods have adopted a basis set of numerical or- 
bitals such as GaussiansSi, LCAG^^'^^i^'' or wavelets^, 
for which the condition of finite support can be exactly 
fulfilled. In contrast, there exists only a few GF transport 
schemes based on the more accurate plane-wave DFT 
methods— This is partly due to the delocalized na- 
ture which excludes the plane waves themselves as basis 
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FIG. 1: Schematic setup used to study phase-coherent elec- 
tron transport. The system is divided into three regions: a 
scattering region (S) and two leads (L) and (R) connecting 
S to thermal reservoirs with chemical potentials fiL and fiR, 
respectively. The electron mean-field potential in the leads 
is periodic and thus all scattering takes place in S. Each 
lead can be build from a principal layer containing an integer 
number of potential periods. 



functions and partly due to the large number of plane 
waves needed to perform even a small calculation with a 
reasonable accuracy. As we demonstrate in this paper, 
both of these problems can be overcome by transforming 
the set of pre-calculated eigenstates into an equivalent 
set of localized Wannier functions. 



C. Conductance formula 

We consider the phase-coherent transport of electrons 
through a system of the form sketched in Fig.^ The sys- 
tem consists of a scattering region (S) connected to ther- 
mal reservoirs via left (L) and right (R) metallic leads. 
The reservoirs have a common temperature but differ- 
ent electro-chemical potentials /il and ^r, respectively. 
We shall assume that the electrons are non-interacting 
and move in a static mean-field potential. The leads 



are assumed to be perfect conductors and thus the elec- 
trons move ballistically in these regions and can scatter 
only on the potential inside S. This assumption has im- 
portant computational consequences which we explore in 
SecHTPl 

By introducing a basis set, {cfii}, consisting of func- 
tions with finite support in the direction of transport we 
can decompose the underlying Hilbert space into disjoint 
subspaces corresponding to the three regions S,L,R. We 
do not require that the basis functions are orthogonal. 
The Hamiltonian matrix, Hij — [(/)_,■), of the entire 

system takes the form 



H 




(1) 



where the Hi themselves are matrices. The vanishing 
coupling between the leads can always be obtained by 
increasing the size of the scattering region, i.e. by in- 
cluding part of the leads in S*. It should be noted that 
if the basis is not orthonormal then in general H will be 
different from the matrix representing the Hamiltonian 
in that basis. 

For non-interacting electrons the retarded single parti- 
cle Green's function can be obtained from the resolvent 
of the Hamiltonian, (zl — H)~^ . We define the Green's 
function matrix by 



(zS-ff)G'-(e) = J, 



(2) 



where Sij — {(f>i\4>j) is the overlap matrix, I the identity 
matrix and z = e + irj, rj being a positive infinitesimal. 
To emphasize the division of the system into the three 
regions we write out Eq. jSJ more explicitly 
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(3) 



It should be noted that G^ differs from the matrix 
{(j)i\{zl — H)^^\(j)j) if the basis is not orthonormal. The 
two are related by 

{cb.\{zI-H)-'\cb,) = [SG"^is)S]^^. (4) 

The Green's function of the scattering region, Gg, plays 
a central role in the theory. It can be writteii^, 

GUe) = {zSs Hs ^l{e) ^'R{e)y' , (5a) 

= (2Ssa-Jfsa)gr(e)(2SL-/fL)(5b) 

= (5c) 

where a € {L,R}. The matrix is a non-hermitian 
self energy which incorporates the coupling to lead a and 



Oa^ie) is the retarded Green's function of the uncoupled 
lead. 

An expression for the current through the system 
was derived by Meir and Wingreea^^ using the non- 
equilibrium Green's function formalism. The most gen- 
eral version of the formula is also valid when interactions 
are present in S. However, here we focus on the non- 
interacting case for which 

2e f 

/=-^ / [nF(e - /ii) - nF(e - mk)] r(e)de, (6) 
where npie) is the Fermi distribution function and the 
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elastic transmission function is given by 

Tie)=Tr[G^sis)TL{e)GUs)TH{e)] 



(7) 



Here Gg — [G^]^ denotes the advanced Green's function 
of the scattering region and Ta is obtained from the self- 
energies as 



r4e) = z(s;(e)-[S;(£)]t) 



(8) 



In the original derivation of Eq. the basis set was 
assumed to be orthonormal, but in fact the formula 
holds for a general basis as has been shown by Xue et 
alP^. Although the current formula jnj is valid for non- 
equilibrium transport we shall focus on the linear re- 
sponse conductance, G. Assuming that the difference 
between and is small and taking the zero temper- 
ature limit we find 

G^GoTiep), (9) 
where the conductance quantum is given by Go = 2e^/h. 



D. Coupling to leads 

The expression Eq. ^ for the self-energy due to lead 
a contains the Green's function g'^^ of the uncoupled 
lead. Since transport through the leads is assumed to be 
ballistic the electron potential must be periodic in these 
regions and consequently the leads can be divided into 
principal layers containing an integer number of potential 
periods. Due to the finite range of the basis functions 
in the transport direction, the size of the principal layers 
can always be chosen so large that only neighboring layers 
couple. In this case the Hamiltonian matrix of the left 
lead takes the form 



H, 



( 



V 



; ; ; \ 

liQ hi 

h\ ho hi 

h\ ho J 



(10) 



where ho is the Hamiltonian matrix of a single principal 
layer and hi is the coupling between neighboring layers. 
A similar matrix structure is obtained for H In writ- 
ing Eq. 1)10(1 we have assumed the lead potential to be 
periodic all the way up to the scattering region. This 
condition can always be fulfilled by extending the scat- 
tering region until it comprises all perturbations arising 
from the presence of scatterers. In practice such pertur- 
bations are rapidly screened by the mobile electrons in 
he metallic leads and the mean field potential is expected 
to decay to its bulk value over a few atomic layers. 

The division of the leads into principal layers with 
nearest neighbor coupling also implies that the scattering 
region only couples to the principal layers immediately 
next to it. From Eq. (jSJ it then follows that only the 



lead Green's function of the first principal layer is needed 
for evaluating the self-energy. This Green's function can 
be determined iteratively due to the periodic structure 
of the matrix {zSa — Ha). A particularly efficient it- 
eration scheme is provided by the so-called decimation 
technique^. 



E. Ultrasoft pseudopotentials 

So far we have made no assumptions regarding the 
form of the mean-field electron potential. In particular it 
could contain non-local terms which are a natural com- 
ponent of most norm-conserving pseudopotentials^SiS. 
However, special care must be taken when ultrasoft pseu- 
dopotentials'^^ are used, as explained below. 

Suppose H is an effective single-particle Hamiltonian 
in which the ion potential has been replaced by an ultra- 
soft pseudopotential. The energy spectrum is found by 
solving the generalized eigenvalue equation 



(11) 



where the pseudo wave functions satisfy the generalized 
orthonormahty condition 



{'4'n\S\ilm) = 5n 



(12) 



Although the Hamiltonian is hermitian its eigenvalues 
do not coincide with the real energy spectrum and it is 
therefore not obvious at the outset how to define the 
Green's function entering the transmission formula (jTji. 

The operator H = S~^H has the correct spectrum but 
it is not hermitian - its eigenvectors are not orthogonal. 
This is, however, not a real problem since self-adjointness 
depends on the inner product. In fact H is hermitian 
if we use the inner product (-I-)! defined by {4'\4>)i — 
{il>\S\(t>). That ('I-)! does define an inner product follows 
from the fact that S is hermitian and positive in the sense 
{(i)\S\(l)) > for all (j). 

The apparent problem related to the use of ultrasoft 
pseudopotentials can thus be removed by changing the 
inner product from (-I-) to (-I-)! and using the Hamil- 
tonian H instead oi H. A direct application of Eq. Q 
leads to the following defining equation for G^: 

{zS-H)-^G''ie) = I, 



where Sij = ((f)i\S\(f)j) and Hij = 



(13) 

\H\(j)j). Conse- 
quently, the only change in the formalism introduced by 
the use of ultrasoft pseudopotentials is the substitution 
of the ordinary overlap matrix {(f)i\(f)j) by the generalized 
overlap {4>i\S\4>j). 



F. Newns Anderson model 

In order to illustrate the use of the general formalism 
introduced above and to develop a physical understand- 
ing of the conductance formula (TJ , we discuss the case of 
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transport through a single energy level coupled to contin- 
uous bands. As we shall see later this model is also useful 
for analysis and interpretation of the transport properties 
of more complicated systems. 

In the Newns Anderson modeSS we consider the sin- 
gle site I a) of energy Sa coupled to infinite leads via the 
matrix elements tan = {a\H\av), where {|q;j')} is a ba- 
sis of lead a. The single level constitutes the scattering 
region and thus all matrix quantities in the transmission 
function Q reduce to complex numbers. The Green's 
function reads 

^-(^) = V (.y (14) 

e- Sa- SL(e) - Sfl,(£) 

with the self-energy due to lead a given by 

Sa(e) = ^W(5a)^.C (15) 

We have dropped the r superscripts since we consider 
only retarded Green's functions in the following. A par- 
ticularly elegant solution to the problem can be obtained 
by introducing the normalized group orbital of lead a: 

\lo.)^V-^Y.^aAav), (16) 

where Vq = (X^iy \^av^Y^'^ ■ The group orbital is located 
on lead a and its weight on a given lead state is deter- 
mined by the magnitude of the corresponding coupling 
matrix element. Consequently, the group orbital is ex- 
pected to be localized at the lead-level interface. The 
coupling between |a) and any state in the lead orthogo- 
nal to 17a) vanishes and thus \a) is coupled to the lead 
via the group orbital only, the coupling being given by 
Va- For simplicity we assume a symmetric contact and 
therefore suppress the a index. Using the general rela- 
tion between a diagonal element of the retarded Green's 
function and the projected density of states (DOS) for 
the corresponding orbital: Im[Gjyi/] = —npi,, we can ex- 
press r (as defined in Eq. Q) by the DOS for the isolated 
lead projected onto the group orbital: F = 7r|Vpp°. Ap- 
plying the same rule to the Green's function of the level 
we arrive at 

r(e) = 27r2|y|2p,(e)p°(£). (17) 

This formula shows that the transmission at a given en- 
ergy depends on three quantities: the coupling strength, 
V , the level DOS, pa, and the DOS of the group orbital 
in the isolated lead, Thus, for an electron with en- 
ergy e to propagate through the system there should be 
states of that energy available in the leads as well on the 
level and the coupling should be reasonably strong. 

An even simpler description is obtained in the so called 
wide band limit where p^J is assumed to be constant. In 
this case the transmission is usually expressed in terms 
of r and Ea and reads 



This is a Lorentzian of unit height and center £„. Its 
width is given by F which also determines the width of 
the energy distribution of \a) in the coupled system and 
thus equals the inverse lifetime of an electron on the level. 
From this expression it is clear that the position of the 
Fermi level relative to Sa is a crucial parameter for the 
conductance. Finally, we note that deviations from the 
simple Lorentzian form H18|l are expected when the cou- 
pling is asymmetric or when p^J (e) has a significant energy 
dependence. 



III. WANNIER FUNCTIONS 

In this section we present a scheme to construct partly 
occupied Wannier Functions (WFs)^. After a brief mo- 
tivation we outline the construction algorithm for both 
isolated and periodic systems. It is natural to introduce 
the algorithm first in the simpler case of an isolated sys- 
tem before extending it to periodic systems, although the 
latter contains the former as a special case. For later use 
we establish an expression for the matrix elements of the 
Hamiltonian in the WF basis. Finally, we illustrate the 
method by constructing the WFs for an isolated ethylene 
molecule and an infinite trans-polyacetylene wire. 



A. Why Wannier functions? 

There are several advantages of using WFs as a basis 
set for transport calculations^: (i) the WFs are spa- 
tially localized which allows for the necessary division 
into scattering region and leads, (ii) any eigenstate be- 
low a certain specified energy can, by construction, be 
exactly reproduced as a linear combination of the WFs 
and thus the accuracy of the original electronic structure 
calculation is retained, (iii) the WF basis set is truly min- 
imal and once constructed the computational cost of the 
subsequent transport calculation is comparable to that 
of semi-empirical methods like tight-binding or extended 
Hiickel. (iv) the WFs contain information about chemical 
properties such as bond types, coordination and electron 
lone pairs and can thus be directly used as an analysis 
tool. 

The WFs are defined as linear combinations of a fixed 
set of single-particle eigenstates - traditionally the occu- 
pied states - with the expansion coefficients optimized 
to give the best localization of the resulting WFs. In 
this paper we shall focus on partly occupied WFs. The 
term "partly occupied" refers to the fact that we allow 
selected unoccupied orbitals in the expansion of the WFs. 
In some situations this can improve the symmetry of the 
resulting orbitals, however, the crucial property that dis- 
tinguish the partly occupied WFs from the traditional 
occupied WFs and make them applicable for transport 
calculations is that they can be localized in metallic sys- 
tems. It is well known that occupied WFs associated 
with a partly filled band in a metal have poor localization 
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propertie3^2i^. The only remedy for this is to include the 
appropriate unoccupied states in the localization space. 



B. Isolated systems 

Consider an isolated system, for which an electronic 
structure calculation produces N eigenstates, {ipn}, of 
which M have an energy below Eo. Our aim is to con- 
struct a set of localized orbitals out of the {V-'n}, which 
span at least the M eigenstates below Eq. We will allow 
eigenstates with energy above Eq in the expansion of the 
WFs in order to improve the localization. 

Before we can start to study localized orbitals, we need 
a measure for the degree of localization. Here we follow 
Marzari and Vanderbilt^ and measure the spread of a 
set of functions, {wn}, by the sum of second moments: 



JV 



{Wn\r\Wn)'^). 



(19) 



We expand the WFs in terms of the M lowest eigenstates 
and L extra degrees of freedom, {(f>i}, from the remaining 
{N — M)-dimensional eigen-subspace. This will lead to a 
total oi M + L WFs. The expansion takes the form 



M 



(20) 



m— 1 



where the extra degrees of freedom (EDF) are written as 



N-M 

E< 

m=l 



(21) 



The columns of the matrix c are taken to be orthonormal 
and represent the coordinates of the EDF with respect 
to the eigenstates with energy above Eq. The matrix 
U is unitary and represents a rotation in the space of 
the functions {-01, . . . , V'Afj 4>It ■ ■ ^ 4'l}- It should be no- 
ticed that the construction H2U|I ensures that the M low- 
est eigenstates are contained in the space spanned by the 
WFs. Through the expansions (|20l) and H21I) . S becomes 
a function of Uij and Cij, which must be minimized un- 
der the constraint that U is unitary and the columns of c 
are orthonormal. The minimization can be performed by 
evaluating the derivatives of 5 with respect to Uij and Cij 
and then using any gradient based optimization scheme 
to iteratively find the minimum. 



C. Ethylene, C2H4 

To illustrate the method we have constructed the WFs 
of an isolated ethylene molecule. We use a cubic super- 
cell of length 14 A and sample the BZ by the F point. 
In Fig. 121 we show iso-surfaces for two different sets of 



WFs corresponding to fully occupied WFs (upper row) 
generated with Eq = Ep and L — Q and partly occupied 
WFs generated with Eq ~ Ep and L = 1. In both cases 
we obtain four equivalent a-orbitals located at the C-H 
bonds. In contrast, the WFs describing the C-C bond 
changes from two equivalent orbitals with a mixed ct/tt 
character, to a single a orbital centered between the C 
atoms and two atomic-like p-orbitals centered on each C. 
By inspection we find that the EDF selected in the min- 
imization procedure in the latter case coincide with the 
anti-bonding molecular 7r-orbital of the C-C bond. Inclu- 
sion of this orbital is precisely what is needed to separate 
the a- and 7r-systems on the molecule. 

Occupied WFs (L=0): 
x4 

x2 




Partly occupied WFs (L=1 
x4 



x2 



© c o o 

00 o ^^^P 



FIG. 2: Occupied WFs (upper row) and partly occupied WFs 
(lower row) for ethylene, C2H4. The symbol xn indicate the 
number of similar WFs obtained at equivalent sites on the 
molecule. Note, that by including the anti-bonding n orbital 
in the WF expansion the o- and 7r-systems become separated. 



D. Periodic systems 

We now turn to the construction of WFs for a periodic 
system. Suppose that {wn} is a set of periodic functions 
defined in a cubic cell of length a. It has been shown 
that in the limit of large a, minimizing the spread (|19|) is 
equivalent to maximizing the functional""' 



IK. 



\Z„ 



where the matrix X is defined as 



X^rn, 



-i(27r/ a)x | 



(22) 



(23) 



The general- 
ization to cells of arbitrar y sh ape is straightforward and 
can be found in Refs. |4lll42l |. Since an isolated system 
can always be studied in a large periodic supercell within 
the F-point approximation, the functional f2 can also be 
used instead of S to construct the WFs of a finite system. 

Since fl corresponds to S only for periodic functions 
in large supercells it seems to be useful only for systems 
that can be treated within the F-point approximation. 
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However, by using the folding properties of the Brillouin 
zone (BZ) this problem can easily be overcome. Consider 
a periodic system with a unit cell defined by the basis vec- 
tors ai , a2 , aa . The union of the k-subspaces associated 
with a uniform Ni x N2 x N3 grid in the l.BZ which 
includes the F-point coincides with the F-subspace of the 
repeated unit cell defined by basis vectors = NiRi. 
Consequently, any function which can be written as a 
sum of functions each characterized by a k from such a 
uniform grid is periodic in the repeated cell and its spread 
can be defined by fl. In this case the matrix elements \2'3\i 
should be evaluated over the repeated cell, however, as 
usual this can be reduced to integrals over the small cell. 
Finally, we note that since the number of k-points deter- 
mines the size of the repeated cell, the original condition 
that the supercell should be sufficiently large is turned 
into a condition of a sufficiently dense k-point sampling. 

As for the non-periodic case, we start from a set of 
single-particle eigenstates, {V'nkjj obtained from a con- 
ventional electronic structure calculation. For a periodic 
system, the eigenstates are Bloch states labeled by a band 
index, n, and a wave vector, k. In accordance with the 
remarks above the k-points should belong to a uniform 
grid which includes the F-point. We denote the total 
number of bands used in the calculation by N, the de- 
sired number of WFs per unit cell by and the number 
of eigenstates at a given k-point with energy below Eq by 
Mk. If Mk > N^, we set Afk — N^. As in the isolated 
case, all eigenstates with energy less than Eq will be con- 
tained in the space spanned by the WFs. We define the 
nth Wannier function belonging to cell i by 



E 



e-^'-^'V^nk, 



(24) 



where Nk is the total number of k-points and V'nk is a 
generalized Bloch state to be defined below. We wish to 
find the set of generalized Bloch states that minimizes 
the spread of the resulting Wannier functions. Since 
Wi^nij) = u'o,ri(r — Rj) it is sufficient to consider only 
the Wannier functions of the cell at the origin. For each 
of the generalized Bloch states we follow closely the idea 
used in the isolated case. Thus we expand ipnu in terms 
of the Mk lowest eigenstates and Lk = — Mk EDF, 
{0/k}; from the remaining space of Bloch states: 

(25) 

m=l ;=i 
where the unoccupied orbitals are written as 

W-A/k 



»;k 



= ^ c(k)„i;'0A/^+„i,k. 



(26) 



m— 1 



Since the states {ipnk} coincide with the F-point eigen- 
states of the repeated unit cell we can use fl to measure 
the spread of Wq,™- The maximization of CI with respect 
to C/(k)y and c(k)y is performed in complete analogy 



with the isolated case. For small systems (< 10 atoms) 
the method is quite stable and usually leads to the same 
global maximum independent of the initial guess used for 
the matrices C/(k) and c(k). In such cases the matrices 
can simply be initialized randomly. For larger systems, 
the method can get stuck in local maxima and we have 
found it useful to start the maximization from an initial 
guess of simple orbitals located either at the atoms or the 
bond centers. 



E. Wannier function Hamiltonian 

The transformation of eigenstates into localized WFs 
induces a transformation of the Hamiltonian from its di- 
agonal spectral representation into a real-space or tight- 
binding like representation. This real-space form of the 
Hamiltonian is in itself very convenient as it provides di- 
rect access to physical parameters for e.g. tight-binding 
models. 

By combining Eqs. H24I25I26|I the nth WF located in 
the unit cell at Ri can be compactly written as 



1 



k,m 



ikRi 



V^(k)„j„V; 



mk 1 



(27) 



where ^(k) can be expressed in terms of J7(k) and c(k). 
The Hamiltonian matrix element connecting the nth WF 
of cell i with the mth WF of cell j then becomes 

iJ(Rj-Ri)„„ = {■W^,n\H\Wj,m) 

= ^e~'(^-^')-''T/*(k),„F(k)i™eik 

(28) 

where ejk = (V'ikl-H^lV'ik) are the eigenvalues. Recall, that 
the WFs as defined here are periodic functions with a 
period given by the repeated cell. This periodicity is re- 
flected in H(R.j — Ri)„,„. In order to describe fully local- 
ized functions the coupling matrix elements must there- 
fore be truncated beyond a cut-off distance given approx- 
imately by Ni/2 unit cells in the direction a^. This means 
in turn that the repeated cell should be large enough, or 
equivalently that the number of k-points should be large 
enough, that the WFs have decayed sufficiently between 
the repeated images. 

To test the accuracy of the transformed WF Hamilto- 
nian we can use it to reproduce the original band struc- 
ture. To this end we form the Bloch basis 



1 



Xnk 



Wi 



(29) 



where Nn is the number of unit cells included in the sum. 
The Hamiltonian matrix, i?(k), in the Bloch basis for a 
given k becomes 



(X„k|i?|Xr„k) =^e*-^^i?(R,)r: 



(30) 



R; 
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The dimension of this matrix equals the number of WFs 
in one unit cell and can be easily diagonalized yielding the 
approximate eigenvalues e„k- The quality of the eigen- 
values depends on the number of cells included in the 
sum (I30|l . i.e. on the range beyond which the coupling 
is truncated. Since iJ(Rj)„m becomes incorrect beyond 
Nil 2 unit cells in direction the quality of the eigen- 
values £„k ultimately depends on the number of k-points 
used in the construction of the WFs. 



F. Trans-Polyacetylene 

As an example of a Wannier function analysis for a 
periodic system we have constructed the WFs of trans- 
polyacetylene and generated the band diagram using the 
corresponding WF Hamiltonian. We use a unit cell con- 
taining two carbon and two hydrogen atoms and with a 
large amount of vacuum separating the repeated wires 
by more than 10 A in the directions perpendicular to 
the wire. We use a uniform 31 x 1 x 1 k-point grid in 
the construction of the WFs. The band diagram of the 
wire obtained directly from the electronic structure cal- 
culation is shown in Fig. 01 (full lines). For transport 
calculations we need a good description of the electronic 
structure of the wire in an energy window around the 
Fermi level, ep — ^- Thus, in addition to the filled bands 
we should describe at least the lower part of the tt* con- 
duction band correctly. Fig. ^ shows iso-surfaces of three 
different WFs obtained using the parameters Nyj = 6 and 
Eq — 1.0 eV. With this choice of parameters we obtain 
six WFs per unit cell and ensure a correct description 
of the band structure up to 1.0 eV. In total we obtain 
a cr-orbital at every C-H and C-C bond and an atomic- 
like p^-orbital centered on each C atom. It is instructive 
to calculate the band structure of the wire using these 
WFs as basis functions as described in section llil El The 
result is shown by dots in Fig. |3| As expected the re- 
constructed bands are in good agreement with the origi- 
nal bands below Eq. The high density of bands starting 
around 3 eV above ep marks the beginning of the con- 
tinuous spectrum. The fact that the reconstructed tt* 
band does not follow a particular band in this region in- 
dicates that the corresponding generalized Bloch states 
(|25|l consist of non-trivial linear combinations of delo- 
calized original bands. This is sometimes referred to as 
band-disentanglement. 
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FIG. 3: Comparison between the band diagram obtained 
directly from the electronic structure calculation (full line) 
and the approximate band structure obtained using the WFs 
of Fig. 0] as basis functions. 
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FIG. 4: Partly occupied WFs for trans-polyacetylene. The 
unit cell contains two CH units. The WFs have been gener- 
ated with the parameters A'^ = 6 and Eo = 1.0 eV, where 
the Fermi level is set to eV. There is one a orbital at every 
C-H and C-C bond and one atomic-like pz orbital on each C. 

of supercells in connection with transport calculations 
we discuss the construction of the three Hamiltonians 

Ha,E[s, Hsa- 



IV. REPRESENTATION OF THE 
HAMILTONIAN 

In this section we link the WFs introduced in sec- 
tion to the general transport formalism of section ITU 
What is needed is a method to construct the Hamilto- 
nian for the coupled L-S-R system in a localized WF 
basis set, and we will discuss a way of achieving this. Af- 
ter a general discussion of some issues related to the use 



A. Periodic supercells 

Within the periodic supercell approach, the system of 
interest is defined inside a finite computational cell, the 
supercell, which is repeated in all directions to form a 
super lattice. The principle is illustrated in Fig. |5l for 
a single molecule suspended between infinitely extended 
surfaces. Clearly, we must require that the transverse 
dimensions of the supercell as well as the thickness of the 
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surface slabs are so large that the periodically repeated 
molecules do not interact. 




FIG. 5: Illustration of a supercell describing a single molecule 
suspended between infinite surfaces. 

Fig. shows an example of a scattering region as de- 
fined in the transport setup sketched in Fig ^ If the 
surface slabs are thick enough, the electron density and 
thus the electron mean field potential at the supercell 
end-planes are close to their bulk values. It is therefore 
permissible to extend the electron potential to the left 
and right of the supercell by the bulk potential which 
leads to the system shown in Fig. We note, that the 
system remains periodic in the directions perpendicular 
to the transport direction. Within this approach, we are 
therefore effectively considering transport through an in- 
finite array of contacts. Whether this provides a good 
description of transport through a single, isolated con- 
tact then depends on the degree of interference between 
the repeated contacts. As the transverse dimensions of 
the supercell are increased this interference is expected to 
decrease and the result should approach that of a single 
contact. 

Since the scattering region of the system shown in 
Fig. El consists of an infinite array of molecules, all ma- 
trices entering the transmission formula have infinite 
dimensions. However, the periodicity of the system im- 
plies that the wave vectors, k^, belonging to the l.BZ. of 
the transverse plane are good quantum numbers. Con- 
sequently, all the matrices entering Eq. Q are diagonal 
with respect to kj^ and the total transmission can thus be 
decomposed into a sum over kj^ -dependent transmissions 
given by 

T(ki ; e) - Tr [G's (k± ; £)rL (kx ; e)G J (k^; e)T R(kr, e)] . 

(31) 

The transmission per supercell is evaluated as T(e) = 
Skj_ ^kj_T(kj_; e), where Wkj_ are appropriate weight 
factors adding up to 1. 

B. Some notation 

For use in the following we introduce some notation. 
We take the direction of transport to coincide with the z- 
axis. For any system of the form illustrated in Figs. I5I()I 
we denote the basis vectors of the super lattice in the 
transverse plane by Ai and A2. We assume that this 
plane is perpendicular to the transport direction, i.e. Ai • 
62 = A2 • =0, and refer to the unit cell spanned 




FIG. 6: The system obtained by extending the electron po- 
tential to the left and right of the supercell shown in Fig.|5]by 
the bulk potential. Note, that the system remains periodic in 
the directions perpendicular to the transport direction. 

by Ai and A2 as the transverse cell. A general lattice 
vector within this plane is denoted by Rj_ , and can thus 
be written Rj_ — nAi -I- mA2, where n,m are integers. 
Finally, we denote the wave- vectors in the l.BZ of the 
transverse cell by kj^. 

C. Lead Hamiltonian {Ha) 

In the following we describe how to construct the 
Hamiltonian matrix, Jfc, for lead a. From Eq. H1U|I it 
is clear that we need the matrices ho and hi describ- 
ing a single principal layer and the coupling between two 
neighboring layers, respectively. 

Suppose that a single principal layer of lead a is de- 
scribed by a supercell with basis vectors A^ = iV^ai, 
i = 1,2,3. If iVi > 1 for some i, the principal layer 
can be build from a smaller unit cell with basis vectors 
ai , a2 , as . In the example illustrated in Fig. [S] a possible 
principal layer super cell consists of one ABC period of 
3x3 sections of the (111) atomic planes of the bulk fee 
crystal. This cell can in turn be made up of smaller cells 
containing only 1x1 sections in the transverse plane. Let 
{wi,n}, n = 1, . . . , N.U] denote the set of WFs obtained for 
the small cell at R^. For a given kj^ wc form the Bloch 
functions 

X^.n{^±) = -7ir= E e'''^'^^^^A^ ' (32) 

where Nji^^ is the number of transverse cells included in 
the sum and i denotes a small unit cell. Note, that even 
though the Bloch functions are extended in the transverse 
directions they preserve the finite range of the WFs in 
the transport direction. The corresponding Hamiltonian 
matrix is given by 

(Xi,n(k_L)|iJ|Xj,m(k±)) 

= J2 e*''^-^^-ff(R, - R. + R±)„™(33) 

where the matrix elements H(Rj — R^ + il±)nm refer to 
the WF Hamiltonian of Eq. if^ . To obtain the matrix 
ho we let i,j run over all small unit cells within one 
principal layer supercell and n, m run over all WFs. To 
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obtain hi we do the same except that j now runs over all 
small unit cells within the supercell of the neighboring 
principal layer. The accuracy of the Hamiltonian H{'k±) 
depends on the number of transverse cells included in the 
sum. As explained in sec. IIII El this number is limited 
by the size of the period of the WFs which in turn is 
determined by the number of k-points used to construct 
them. For a transverse cell containing 3x3 atoms we have 
found it sufficient to sum over nearest neighbor cells. 



D. Hamiltonian for scattering region (Hs) and 
coupling (Hsa) 

The technique for constructing a matrix representation 
of the lead Hamiltonian described in the previous section 
can also be applied to obtain the Hamiltonian for the 
scattering region, Hs- However, the coupling between S 
and the leads, Hsa, involves matrix elements connect- 
ing WFs in the lead with WFs in S, and since these are 
in general different the method described above does not 
apply to the coupling matrix. Instead, we use the elec- 
tronic structure code to directly evaluate these matrix 
elements. 

We define a composite supercell, S, by attaching one 
lead principal layer cell on each side of the supercell con- 
taining S. The electron density and the atomic config- 
uration in the new supercell are obtained by "gluing" 
together the densities and atomic configurations of the 
three constituent parts. Since the electron density defines 
the effective potential, this determines the Hamiltonian 

Let us denote by Bi the combined set of WFs for S 
and the two principal layers closest to S. For each WF 
in Bi we form the Bloch function (|32|l on a real space grid 
in the composite supercell S. The resulting set of Bloch 
fmictions is denoted B2. Note, that the Bloch functions 
in the lead parts of S are identical to those used to con- 
struct the lead Hamiltonian. Since the WFs are originally 
defined in separate simulation cells, the introduction of 
these WFs into the composite cell S involve some manip- 
ulations which are easiest to implement in real space. For 
illustration, consider a WF of S constructed using only 
a single k-point in the transport direction. Since the pe- 
riod of the WF is given by the length of S, it cannot 
be directly introduced in S. Instead, we first translate 
it to the center of S (using periodic boundary conditions 
within S) , then copy it into the S part of S and extend it 
into the lead parts of S by zero-" tails" . Finally, the WF 
is translated back to its original position (using periodic 
boundary conditions within S). 

Given the set of Bloch functions, B2, defined in the 
composite supercell as well as the Hamiltonian, Hg, we 
use the electronic structure code to evaluate the Hamil- 
tonian matrix elements. This results in a matrix of the 



form 

/ ho Hl^ hi \ 
Hg = HsL Hs HsR ■ (34) 

V hi hIj, ho J 

The matrices hi and h[ are due to the periodic boundary 
conditions. Except for a constant, the matrix ho should 
equal ho and a comparison of the two can serve as a con- 
sistency check. Finally, Hg must be shifted by a constant 
to account for the unspecified energy-zero and thereby 
ensure a smooth matching at the interface between lead 
and scattering region. This can be done by adding to 
Hg the matrix cSg, where c = [/i.o]oo — [^o]oo- 

V. RESULTS AND ANALYSIS 

In this section we apply the transport scheme to two 
molecular contacts: a hydrogen molecule between atomic 
Ft chains and benzene-dithiol between Au(lll) surfaces. 
We introduce a useful technique for relating the various 
features of the transmission function to specific molecular 
orbitals. This enables us to identify the anti-bonding H2 
orbital as the current carrying state, and to explain the 
transmission through the benzene-dithiol in terms of only 
two molecular orbitals. 



A. Molecular hydrogen bridge in a Pt chain 

As a first example we study the transmission through 
a molecular hydrogen bridge between a pair of semi- 
infinite, monatomic Pt chains, see Fig. Ela). In a re- 
cent experiment the conductance of a molecular hydro- 
gen bridge suspended between bulk Pt electrodes was 
measured using mechanically controlled break junctionsi. 
The system was found to have a conductance close to the 
conductance quantum. Go — 2e'^ /h. While this has been 
confirmed by at least three different calculationsii^S^i, 
there is at least one calculation predicting a conductance 
of only 0.2Gg^^. For simplicity we use monatomic chain 
leads where the transverse k-point sampling can be safely 
omitted. 

We represent the system in a supercell with trans- 
verse dimensions lOAxloA. The relevant bond lengths 
are: dpt-pt = 2.41 A, dpt-H = 1-67 A and dn-H = 0.88 A. 
We include 4 Pt atoms in a principal layer and define the 
scattering region as 4Pt— H2— 4Pt. 

In Fig. |H1 we show the calculated transmission (thick 
line) as a function of energy. The transmission vanishes 
below —6.5 eV which marks the bottom of the s-band of 
the infinite Pt chain. The predicted conductance, G — 
T{eF), is found to be close to the experimental value of 
IGq. Ideally, we would like to understand the shape of 
the transmission curve in terms of the properties of the 
uncoupled systems, i.e. the isolated H2 molecule and the 
free leads. The interaction of a molecule with the leads 
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(a) 





(b) 



f H2 anti-bonding state f 

Pt d-orbital / . , 

Pt o -orbital 

FIG. 7; (a) Molecular hydrogen bridge connecting two semi- 
infinite, monatomic Pt chains. The atoms shown constitute 
the scattering region in the conductance calculation, (b) Con- 
tour plots of three representative basis functions. After diag- 
onalizing the molecular subspace, the basis set consists of the 
initial WFs in the lead region and renormalized H2 bonding 
and anti-bonding states in the molecular region. 



can be seen as a two-step process: first the molecular 
levels are shifted (renormalized) due to the change in the 
potential from the leads. Next, the renormalized levels 
hybridize with the states in the leads as a consequence of 
the overlap between the wave functions. In the following 
we describe how the first step can be implemented and 
studied in practice. We will return to the second step in 
Sec.lVRl 



1 


4- 


1 


2- 









1 


en 




|o 


8- 








6 




4- 





2- 




0- 



Total transmission 

bonding state removed 

Hg anti-bonding state removed 




-4 -2 
Energy (eV) 



FIG. 8: Transmission through the molecular hydrogen bridge 
shown in Fig. |7[a). The thick line represents the total trans- 
mission while the dashed and thin lines refer to the situations 
where the H2 bonding, respectively anti-bonding states have 
been removed. 



1. Diagonalization of the molecular subspace 

In order to obtain the renormalized energy levels and 
orbitals of a molecule in contact with leads, we first con- 



struct the WFs for the combined system, i.e. the scat- 
tering region, and the corresponding Hamiltonian ma- 
trix, Hs- Next, we identify the set of WFs located on 
the molecule. This involves some arbitrariness since it 
is not always clear where to put the separation between 
molecular WFs and lead WFs. Usually, the distinction 
is clearer when the bonding between molecule and lead 
is weaker. We refer to the space spanned by the molec- 
ular WFs as the molecular subspace. The renormalized 
molecular levels and orbitals are obtained by diagonaliz- 
ing Hg within the molecular subspace. Since molecular 
orbitals to a large extent are determined by the symmetry 
of the molecule, it is usually straightforward to relate the 
renormalized orbitals to the orbitals of the free molecule. 
In contrast, the renormalized energies can be drastically 
shifted compared to the energies of the free molecule. 

In the case of H2 the molecular diagonalization pro- 
duces a bonding, and anti-bonding, |a), state from 
the initial Is-like WFs centered on each H. The renor- 
malized energies are Sb — {b\H\b) = —7.4 eV and Sa = 
{a\H\a) = 0.4 eV, relative to the Fermi level of the in- 
finite chains. For comparison the corresponding energy 
levels of the free molecule before coupling are —4.2 eV 
and 6.6 eV, respectively. The anti-bonding state is thus 
shifted down close to the Fermi level. This agrees with 
the conventional understanding of hydrogen dissociation 
on transition metal surfaces which has been established 
on the basis of DFT calculations''^. In fact the filling 
of the H2 anti-bonding state weakens the hydrogen bond 
and eventually leads to dissociation. The effect is less 
dramatic for hydrogen in the contact, however, the frac- 
tional filling of the anti-bonding state does enlarge the 
H-H bond length as compared to the free molecule. The 
reduction of the HOMO-LUMO gap from 10.8 eV in the 
free molecule to 7.8 eV in the contacted molecule is partly 
due to this structural effect (~ 1 eV). The remaining re- 
duction is due to the lead-induced change of the potential 
around the molecule. In Fig.jjfb) we show contour plots 
of \a) together with Pt d- and cr-like WFs. 

Having obtained the renormalized molecular orbitals, 
it is natural to ask if there is significant interference be- 
tween the two levels or if the conduction properties of the 
system are determined by just one of the two orbitals. A 
very direct way to answer this question, is to calculate 
the transmission with cither the bonding or anti-bonding 
state removed. In practice this is done by setting all 
Hamiltonian matrix elements involving either \b) or \a) 
to zero. The resulting transmission curves are shown in 
Fig. |S1 It is clear that the transmission is mainly due to 
the anti-bonding state and that interference effects are 
not very pronounced. The small peak in the transmis- 
sion through the bonding state around the Fermi level is 
due to a peak at the same position in the DOS of the 
group orbital of \b) (see Sec. IIIF|l . Although the use of 
monatomic chains as leads is an over-simplification, the 
main conclusions regarding the value of the conductance, 
the position of the molecular levels and the mechanism of 
conduction are unchanged when realistic bulk electrodes 
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are usedai. 



B. Benzene-dithiol between Au(lll) surfaces 

In this section we study the transmission through 
a benzene-dithiol molecule (BDT) suspended between 
Au(lll) surfaces. The transport properties of this sys- 
tem have been investigated experimentally by Reed et 
al.^ using the mechanically controllable break junction 
technique, and several theoretical studies have been re- 
ported subsequentljii^*22iiLi2iiS. All of these studies are, 
like the present, based on a combination of the Landauer 
Biittiker formalism where a DFT description of the elec- 
tronic structure enters at different levels. In general, the 
calculations for a BDT molecule covalently bonded to 
the electrodes finds a conductance which is 2-3 orders 
of magnitude higher than the experimental value. Dif- 
ferent reasons for this discrepancy have been proposed 
including overlapping BDT molecules^ and attachment 
of the BDT on top of a gold atomic. Here we focus on 
a symmetrically, covalently bonded molecule. We find 
a conductance comparable to those obtained in previous 
theoretical studies, however, we illustrate in addition how 
a proper construction of renormalized molecular orbitals 
leads to a very simple picture of the transport mecha- 
nism in terms of only two molecular states. The rela- 
tively large number of conductance calculations reported 
for the BDT system makes it an ideal test system for 
comparing different computation methods. 

We use a supercell containing 3x3 Au atoms in the 
directions perpendicular to the transport. The lead prin- 
cipal layer contains three atomic planes in accordance 
with the ABC-stacking. The scattering region contains 
3 Au layers on each side of the DTB molecule which is 
enough to achieve a smooth matching with the bulk po- 
tential. To obtain a contact geometry which more closely 
resembles the open structures likely to occur in the ex- 
periments we insert a 3-atom pyramid between the sur- 
face and the molecule. The central part of the scattering 
region is shown in Fig. IHta) . We use the bond lengths 
dAu-S = 2.40 A, ds-c = 2.75 A, dc-c = 1-40 A and 
dc-H — 1-08 A in accordance with those used in Ref. 14^. 
All gold atoms are fixed in an fee lattice with lattice pa- 
rameter a = 4.18 A. 

In Fig. ^1 we show the calculated transmission ob- 
tained using 1, 8 and 18 irreducible k-points to sam- 
ple the transverse BZ. The F-point transmission fluctu- 
ates very much compared to the converged result and 
is clearly not a good approximation. The rapid fluctua- 
tions are characteristic of F-point transmissions. This is 
because the lead in this case is simply a one-dimensional 
pipe with a 3 X 3-atom cross section. The electronic struc- 
ture of such relatively thin Id structures is very different 
from the electronic structure of bulk systems, in partic- 
ular the DOS of the former exhibits van Hove singulari- 
ties at the Id band edges^*'. We mention that the band 
dispersion between the free BDT molecules is less than 
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FIG. 9: (a) Benzene-dithiol (BDT) between Au(lll) sur- 
faces, (b) The three renormalized molecular orbitals closest 
to the Fermi level of the electrodes. The renormalized energy 
is indicated in parenthesis. Orbital 2 (HOMO) has negligible 
weight on the S atoms and has no influence on the conduction 
properties. In contrast orbitals 1 and 3 are strongly coupled 
to the lead states and in fact they completely determine the 
transmission function. 



0.1 eV for all relevant orbitals. This shows that the k- 
point dependence of the transmission function is not due 
to direct coupling between the repeated molecules, but is 
rather a substrate induced effect. The calculated trans- 
mission shows two broad peaks at around —3.0 eV and 
2.5 eV, and two more narrow peaks around —1.5 eV and 
— 1.0 eV (all energies are with respect to the Fermi level 
of the leads). The occurrence of two broad peaks, one 
below and one above the Fermi level, is in qualitative 
agreement with previous results^ii^S*^. The position of 
these peaks as well as the existence of the two narrow 
peaks just below the Fermi level, are not in agreement 
with the earlier calculations which also deviate somewhat 
from each other. These discrepancies might be due to 
differences in the adopted geometries or technical issues 
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like ademiate k-point sampling. In particular the result 
in Ref.|43has been averaged over different atomic geome- 
tries which might smear out features related to details in 
the lead electronic structure. As we shall see, the narrow 
peaks arise exactly due to such details. Finally we note, 
that the predicted conductance is around 0.1 Go- This 
is in very good agreement with the results obtained in 
Refs. 47 49^ while it is a factor of 4 smaller than the re- 
sult of Ref. li^ and a factor of 3 larger than the result of 
Ref.il^ 



clearly demonstrates that the broad peak at —3 eV is 
exclusively due to orbital 1, while the peak at 2.5 eV is 
exclusively due to orbital 3. In particular interference 
between the different molecular orbitals does not occur 
around these energies. We can also conclude that the 
narrow peak at —1.5 eV is mostly due to orbital 3 while 
the narrow peak at —1.0 eV is mostly due to orbital 1, 
although the separation is not so clear in this case and 
interference effects do play a role. 
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FIG. 10: Calculated transmission through the bdt-system 
shown in Fig. El for three different k-point samplings of the 
transverse BZ. 

To investigate which orbitals of the BDT molecule 
are responsible for the different features of the transmis- 
sion function we have performed a diagonalization of the 
molecular subspace as described in Sec. lV A Tl We define 
the molecular subspace as those WFs whose centers lie 
closer to the C6H4 molecule than to any other atom. The 
reason why we do not include the WFs at the S atoms 
in the molecular subspace is that S is three-fold coordi- 
nated to the gold and the "weak" link in the contact is 
more naturally defined between S and C. Diagonalizing 
the Hamiltonian within the molecular subspace leads to 
16 renormalized orbitals distributed in the energy range 
-17 eV to 5 eV. In Fig.|^b) we show three renormalized 
molecular orbitals with energy close to the Fermi level. 
These orbitals all belong to the 7r-system of the molecule, 
i.e. they change sign under a reflection in the molecular 
plane. 

Orbital 2 has negligible weight in the vicinity of the 
S atom and is expected to couple only weakly to the 
lead states. In contrast orbitals 1 and 3 have signifi- 
cant weight at the S-C bond and should therefore couple 
strongly with the leads and contribute correspondingly to 
the transmission. To obtain more quantitative informa- 
tion we have calculated the transmission function with 
all molecular orbitals except orbital 1, respectively 3, re- 
moved from the basis set, see Fig. The result very 
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FIG. 11: Calculated transmission through the bdt-system 
when all molecular orbitals except from orbital 1 (thin) and 
orbital 3 (dashed) have been removed. The total transmission 
(thick) is shown as a reference. All calculations have been 
performed using 8 irreducible k-points. Notice, that the main 
transmission peak at —3 eV (2 eV) can be directly related to 
orbital 1 (3). 

The above analysis shows that the transport properties 
of the BDT are completely determined by orbitals 1 and 
3. To take the analysis even further we have constructed 
the corresponding group orbitals (see Sec. Ill F|) . both of 
which are quite similar to an S-centered p-orbital perpen- 
dicular to the molecular plane. This is easily anticipated 
from the tt character of the molecular orbitals. Accord- 
ing to Sec. Ill Fl the DOS of the group orbital calculated 
without coupling to the molecule, together with the level 
position and the coupling matrix element, determines the 
projected DOS for the corresponding molecular orbital. 
The group orbital DOS is shown in the upper panel of 
Fig. El (dashed hne). We can model the DOS by a semi- 
elliptical band on top of a flat background (solid line), 
which allows us to compute the resulting DOS of an or- 
bital with energy Sa and coupling V. The coupling can 
be directly extracted from the DFT Hamiltonian giving 
V = 1.3 eV. Using the on-site energies of the renormal- 
ized orbitals 3 and 1, we obtain the DOS shown in the 
middle and bottom panel, respectively. In addition to 
the broad resonances formed close to the bare on-site en- 
ergies, a narrow bonding, respectively anti-bonding reso- 
nance, involving the molecular orbital and the lead states, 
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is formed at the edges of the semi-elhptical band. On the 
basis of these observations and the conductance formula 
Eq. H17|l we can finally conclude that the transmission 
peaks at 2.5 eV and —1.5 eV are due to orbital 3, while 
the peaks at —3.0 eV and —1.0 eV are due to orbital 1. 

It is worth noting that the simplicity of the picture 
described above in terms of only two molecular orbitals 
relies on ou r definition of the molecular subspace. In 
Refs. I29l48l the transmission was analyzed in terms of 
the molecular orbitals of the whole BDT molecule^, i.e. 
including the S atoms, and the BDT molecule plus se- 
lected Au atoms^^. This leads to a more complicated 
description of the transmission peaks involving several 
different molecular orbitals. Thus the precise division 
of the Hilbert space into molecule- and lead-subspaces 
can lead to more or less elegant descriptions of transport 
properties. 



periodic systems and examples have been given to illus- 
trate the important feature of bonding-antibonding clo- 
sure. Finally, we have applied the transport scheme to a 
molecular hydrogen contact in a monatomic Pt wire and 
a benzene-dithiol molecule between Au(lll) surfaces. A 
useful analysis technique for identifying which molecu- 
lar orbitals contribute to which features of the transmis- 
sion function has been introduced and applied to each of 
the systems. In this way we showed that the transport 
through the hydrogen molecule is determined by the anti- 
bonding state, and we identified two molecular orbitals 
which completely accounts for the transport through the 
benzene-dithiol. 



VII. ACKNOWLEDGMENTS 

We thank Lars B. Hansen for contributions to the Wan- 
nier function scheme. We acknowledge support from the 
Danish Center for Scientific Computing through Grant 
No. HDW-1 101-05. 




FIG. 12: The upper panel shows the DOS, p2^, of the (com- 
mon) group orbital of molecular orbitals 1 and 3. has 
been calculated with all coupling matrix elements involving 
the molecular orbitals set to 0, i.e. it refers to the uncoupled 
leads. The dashed line is the result obtained from the calcu- 
lation, while the full line is a model fit. The two lower panels 
show the resulting DOS of a level 



VI. CONCLUSIONS 



We have presented a numerical method for calculat- 
ing electron transport in atomic-scale contacts. The 
method combines a formally exact Green's function for- 
malism with a mean-field description of the electronic 
structure based the Kohn-Sham scheme of density func- 
tional theory. By transforming the delocalized eigen- 
states obtained from an accurate plane-wave calculation, 
into maximally localized Wannier functions we obtain a 
highly efficient minimal basis set for evaluating the rele- 
vant Green's functions. The construction scheme used to 
obtain the WFs has been introduced for both isolated and 
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